Spin Nematic Phase in S—l Triangular Antiferromagnets 
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Spin nematic order is investigated for a S—l spin model on triangular lattice with bilinear- 
biquadratic interactions. We particularly studied an antiferro nematic order phase with three- 
sublattice structure, and magnetic properties are calculated at zero temperature by means of 
bosonization. Two types of bosonic excitations are found. One is a gapless excitation with lin- 
ear energy dispersion around k ~ 0, and this leads to a finite spin susceptibility at T = and 
would have a specific heat C(T) ~ T 2 at low temperatures. These behaviors can explain many of 
characteristic features of recently discovered spin liquid state in the triangular magnet, NiGa2S4. 



The concept of spin liquid was introduced thirty years 
ago by P. W. Anderson |lj, as a quantum critical state 
in which the spin-spin correlation function does not show 
a real long-range order but a power-law behavior. This 
issue has been studied intensively since then both theo- 
retically and experimentally. Frustration and quantum 
fluctuations are considered two ingredients to realize a 
spin liquid, and a spin-1/2 Heisenbcrg antiferromagnet 
on triangular lattice was the first candidate. Many anti- 
ferromagnetic materials with triangular lattice structure 
have been studied to see if they may show a spin liq- 
uid behavior, but most of them turned out to exhibit 
some long-range order at low temperatures. Very few 
exceptions are 3 He thin layer Q and organic k-(BEDT- 
TTF) 2 Cu2(CN) 3 0, and just recently a new spin- liquid 
material NiGa2S4 was discovered While the former 
two are spin-1/2 system, NiGa2S4 is a spin-1 system. 

In NiGa 2 S4, spins of Ni 2+ ions form triangular layers 
with undistorted regular triangle units, and the layers 
are stacked along the c-axis. These layers are effectively 
decoupled, since the Ni-Ni distance is more than three 
times longer between layers. This system showed vari- 
ous low-temperature properties that indicate a spin liquid 
state. First of all, no singularity was observed in specific 
heat down to the lowest temperature, T—0.3 K, meaning 
the absence of phase transitions. Moreover, the specific 
heat shows a power-law behavior, C ~ T 2 , below 10K. 
Secondly, the magnetic susceptibility gradually increased 
with decreasing temperature, and approached a finite 
value. Thirdly, neutron experiment observed a peak at an 
incommensurate wavevector Q~(-2=,0). However, this 
was not a magnetic Bragg peak, and spin correlation 
length did not diverge but saturated to about £~ 20A, 
only seven lattice units. 

Absence of magnetic long-range order and presence of 
critical behaviors are necessary conditions to identify a 
spin liquid, and these were satisfied in NiGa2S4. These 
may suggest a finite spin gap instead, but this contradicts 
a nonvanishing temperature dependence of susceptibility. 
In this paper, we will examine the possibility of a hidden 
order that reproduce similar low-temperature properties 
as critical spin liquid states. 

Possible order parameters are not ordinary static spin 



FIG. 1: Spin quadrupole moment Q M/J ' of a single-site wave- 
function when (a) (S) / and (b) (S) = 0. (c) Three- 
sublattice nematic order. Dotted triangle shows a unit cell of 
the ordered state. 



dipole moments (S), since neutron experiment did not 
observe magnetic Bragg peaks. We should note that 
Ni 2+ ions do not have orbital degrees of freedom p|, 
and that we can describe this system as a pure spin 
model with no spin anisotropy. Therefore, if any long- 
range order exists, its order parameter should be rep- 
resented in terms of spin operators. We will investi- 
gate the simplest candidate, spin quadrupole moments, 
Qw*' = |(^^V^'^)-i5(5+l)5 w /, where /x is spin 
index, and this also corresponds to nematic order 0,|(J. 
The nematic order parameter Qn V describes anisotropy of 
spin fluctuations, not static moment, and can be nonzero 
only if S > 1 0). In NiGa2S4, local spins are S=l and 
therefore we consider the nematic order parameter de- 
fined at each site. Neutron experiment found a peak of 
scattering at incommensurate wavevectors Q, not at q=0 
or Brillouin zone boundary. This suggests that the ex- 
pected nematic order is not uniform but modulates in 
space, i.e., some antiferro order. 

To describe this order, the standard Heisenberg Hamil- 
tonian is not sufficient, since it is believed to have a 120- 
degree magnetic long range order |||. We use a spin-1 
model with additional biquadratic interactions between 
nearest neighbor sites on the triangular lattice, 



H = 



(i,j) 



[JSi-Sj+KiSi-Sjf] 
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This model has been studied by mean field analysis and 
numerical calculations 0- The mean field analysis at 
T = showed two nematic phases. One is the region 
of -KT<J<0 and a uniform nematic order appears. The 
other one is an antiferro nematic order and it is predicted 
in the region of 0< J<K. The latter is our expected case, 
and we will investigate that parameter region. We should 
note that the model (JTJ is a phenomenological Hamilto- 
nian introduced to describe an expected modulated ne- 
matic order. While the antiferromagnetic bilinear terms 
are naturally expected from supercxchange processes, the 
biquadratic terms are also present as higher-oder pro- 
cesses of virtual electron hopping, but usually it is ex- 
pected that |J|>|if|. In our choice of model, it is as- 
sumed that both coupling constants are renormalized in 
the low-energy sector from their microscopic values such 
that 0<J<K, but this should be examined in future 
work. 

We investigate the nature of nematic order parame- 
ter <3 p „ in more details for the S=l case. To see this, 
it is more convenient not to subtract the constant term 
in the definition, and consider Q lll/ (r) = ^(S >J '(r)S 1 ' (r) 
+S v (r)S >J '(r)) at site r. We introduce this tensor at 
each site r defined for a local wavefunction, and calcu- 
late its average by taking account of fluctuations in space 
and time. In the 5 = 1 case, we can show that for any 
single-site wave functions, the three eigenvalues of Q^ u 
are Ai,2=|(l T \A (S) 2 ), and A3=l, and Q^ v is repre- 
sented as an ellipsoid as shown in Fig. QJa). Therefore, 
when (S)=0, spin fluctuations are like disk and have zero 
amplitude in one direction in spin space (see Fig.^b)). 
This is characterized by the director vector n that is 
perpendicular to the disk. The corresponding wavefunc- 
tion |i/) n ) is the coherence state such that (n • S)|^ n )=0, 
namely rotated \S Z =0) state with quantization axis par- 
allel to n. 

Mean field solution of the model QJ is easily obtained 
for the triangle lattice, and it is a three-sublattice ne- 
matic order shown in Fig.^c). The directors in the three 
sublattices are orthogonal to each other, and without 
generality, we can set n^||a;, ns||j/, and nc||z, where A-C 
are sublattice indices. Then the mean-field ground state 
is a direct product of local states | x I / mf)= IIr \S x —0)a,r 




FIG. 2: Two types of excitations for J/K=0. 5. (a) Density 
of states with logarithmically divergent van Hove singularity, 
(b) Energy dispersion. 



S CR- a CR+ a CR, S CR-^CR+^CR. S C R --z(a:] 7R /3cR 

—(3^ CR acR)- In a similar way, two types of bosons are 
also introduced for the A- and B-sublattices. It is noted 
that these bosons are subject to constraint, a] R OjR 

+PjR0jR< 1 > or equivalently a 2 R =/3 2 R =a jR (3j R =0. 

We replace spin operators in the model JIJ by bosons 
and then neglect their interactions. This corresponds 
to Gaussian approximation of spin fluctuations. After 
taking the Fourier transformation, the obtained boson 
hamiltonian reads, Hb=3KQ+ H{(3a,chb)+ H{(3b,(xc)-\- 
H{(3 c ,oca), and 

H((3j, oyO = 3K J2 (/*Jkftk + «k^' k 

k 

+ 3 E [(J-K)0} k a j/k + J/3] k a], k ] + h.c.} (2) 

k 

where f2 is the number of unit cells and the sum is 
taken over the reduced Brilloum zone of the three- 
sublattice order, and 7k=^e _4fc;c + §e lfc */ 2 cos(v / 3fc ?/ /2). 
Here we set the lattice constant being unity. Each 
type of bosons interact with only another one type of 
bosons, and therefore we can easily diagonalize the bo- 
son hamiltonian by Bogoliubov transformation, Hb = 



ergy 



E 



m=± ' 



k i ra k + Eq with the eigenen- 



\Sy=0)B,R 



S z =0) C ,R. Here, r=(j, R), with sublat- = 3X V (1 ± W)(l ± «fok|), « = 1 - 2J/#. (3) 



tice index j € {A, B, C} and unit-cell coordinate R. 

We now describe quantum fluctuations of the mean- 
field solution by introducing bosons. Each S—l spin 
has three local states, and one of them corresponds to a 
mean- field solution. Following Matveev jfj, we consider 
it as a local vacuum and the other two excited states as 
one-boson states. To be specific, let us consider a site 
in the C-sublattice. \S z =0)c,r is local vacuum |vac)c,R, 
and the other states are represented as \S z =±1)c,r= 
2 -1 / 2 (a^ R ±?/3^, R ) |vac)c,R, by boson creation opera- 
tors cr and $ . Then spin operators are represented as 



The dispersions and density of states are plotted in 
Fig. [21 for J/K =0.5. The £_ branch is gapless ex- 
citation with asymptotically linear dispersion, £_.k ~ 
3-v/ JK/2 |k| around k=0. The e + branch is gapful ex- 
citation, and it touches the £_ branch at energy 3K on 
the six corners of the Brillouin zone. The J = case 
is special, since the mean-field state is an exact eigen- 
state, and the gapless excitations have a quadratic dis- 
persion, e_k oc k 2 . The J—K case is also special and 
the e+ branch becomes gapless and degenerate with the 
e_ branch. 
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FIG. 3: Dynamical spin structure factor y\ S^(k, ui) for 
J/K=0.5. Dispersion of two magnon branches is also plotted. 
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FIG. 4: Static structure factor S(k) for J/K=0.5. 
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These two branches describe bosonic elementary ex- 
citations in the nematic order, and in particular, the 
gapless branch is the Goldstone mode corresponding to 
the broken spin rotation symmetry [Io| . These elemen- 
tary excitations contribute to magnetic fluctuations, and 
therefore, we may call them magnons also in this case. 
The Gaussian approximation is checked by calculating 
p= (ala r +(3l/3 r } , the density of mean-field excited states. 
p(J) increases from at J=0 to about 0.15 at J=K, 
which was much smaller than 1, and the approximation 
is justified semi-quantitatively. 

The dynamical spin structure factor is given at zero 
temperature by sg (k,u,)= E,<0|^_» H^' k |0) 
5(lu— Ev+Eq), where \v) is the eigenstate with energy 
E v and |0) is the ground state. Sj k is the Fourier trans- 
form of the spin on the j-sublattice. Structure factor 
SjP (k, uj) is diagonal with respect to spin indices p and 
pit. {S^,{k,oj)} is obtained as 
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D AA D AB 
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D BB 
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(0) 





-i<t>k pv 1 ) 



(o) 



(4) 

where e^ k = 7k/|7k|- ^(k,^) and F 2 (k,uj) denote 
one and two-magnon contributions: 



F 2 (k, 



\ ^(-mfe 2e -^- £m , k ), (5) 

m— ± 

sinh 2 (6» m , k+q - 0„, q ) 

q,m,n— ± 
X5(UJ — £ m ,k+q — £n,q), 



0.4 - 



0.2 - 



.total 



2-magnon (x20) 
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FIG. 5: Magnetic susceptibility \m- Two-magnon contribu- 
tion is multiplied by factor 20. 



each k, there are two delta-function peaks at the one- 
magnon energies £±,k, and they are accompanied by two- 
magnon continuum on the higher energy side. As k — > 0, 
the delta- function peak of the gapless branch vanishes as 
\jKj (8J)|k|<S(tj — £_,k)- There are no magnetic Bragg 
peaks, consistent with the absence of ordinary magnetic 
dipole order. 

In Fig. 0] we present the static structure factor S(k) = 
E r (0|5^(r)S'' i (O)|0)e lk r for J/K = 0.5. In the the orig- 
inal Brillouin zone of the triangular Bravais lattice, we 
have 



S(k) = | ^ (1 — TOCOStftkJe" 

m— ± 

+ 1272 Sinh2 ( m,k+q 
q,m,n— ± 



(7) 



Note that we have </>k+t> ± = </>k + 27r/3, where b± = 
(^ji- 2 ^). The one-magnon contribution of the gapful 
mode is strongly suppressed near k = 0, since we have 
cos^k ~ 1 - k 2 x (kl - 3fc 2 ) 2 /1152. Around the iG-point 
(6) k = |(27r,27r/\/3), we have e 2e ±- k - 1 - J/(2K)\k 



k |, and the S(k) shows a singularity of cone tip. This 
leads to a power-law decay in the real-space spin-spin 
correlation as |r| -3 . 

The magnetic susceptibility x-m is calculated from 
the dynamical correlation function using the re- 
lation, Xm=lim k ^o §E M 'e{A,B,c} Jo° duS$(k,u})/w. 
tively. Figure 01 shows ^ - S^(k, u>) for J/K = 0.5. For As k— >0, the one-magnon contribution of the gapful 



where e 4m6 'm,k— (\ _|_ mK|7k|)/(l + w|7k|)- The summa- 
tion of the momentum q is taken over the reduced Bril- 
louin zone in the two-magnon contribution. {S'JJ, (k, lo)} 
and {Sjp (k, w)} are given by replacing sublattice indices 
(A,B,C) in Eq.0) by (B, C, A) and (C,A,B), respec- 
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mode vanishes but that of the gapless mode converges 
to a finite value, 2/(9 J). Adding the two-magnon contri- 
bution, the result is 



2 

9J 



1 ^sinh 2 (0 +jq -0-, q ) 



q 



£ +,i ' 



(8) 



in physical units (g/is) 2 . It is noted that Xm is isotropic 
in spin space. The one-magnon part is independent of 
K, and agrees with the classical value calculated by the 
mean-field approximation. This value also coincides with 
the classical value for the 120-degree order in the pure 
Heisenberg model. In Fig. |SJ J-dependence of Xm is 
shown. It is dominated by the one-magnon part and the 
two-magnon part is very small, about 1.85% at largest. In 
both J=0 and J—K cases, the two-magnon part of Xm 
vanishes, and ex J around J —0 and oc\/K — J around 
J = K. 

Similarly we can calculate the nematic correlation 
(Q^( r )2/-«'f'( u ))- This tensor can be decomposed 
into (1 - <S M „)(1 - V"')(W^"' + <W*«f*')^$(r) + 



, (r). The Fourier transform G^ (k) has a 
similar structure as in Eq.Q, which consists of one 
and two-magnon parts, while G K (k) has the delta- 
function peaks reflecting the existence of static nematic 
quadrupole moments and two-magnon part. The one- 
magnon part of the gapless mode in G$ (k) diverges as 
|k| _1 in the limit k — > 0. The two-magnon part diverges 
more slowly as log |k| around k = in both G^i;(k) and 
G^(k). 

Let us discuss the implications of the above results 
and compare with the experiments for NiGa2S4 0|. We 
have found gapless bosonic excitations with linear en- 
ergy dispersion. They contribute to specific heat as, 
G~ 127rC(3) (T/t;) 2 ~ 45.3(T/v) 2 in units of k B per spin, 
where v—y / 9JK/2 is the velocity of the gapless mode. 
Since the order parameter is a tensor with continuos de- 
grees of freedom as for the 120-degree magnetic order 
[ll|, the nematic order does not appear at finite tem- 
peratures in two-dimensional systems [l^ . The zero- 
temperature magnetic susceptibility is finite and given 
by Xm ~ 2/ (9 J). These behaviors agree with the ex- 
perimental data. The spin structure factor does not 
show magnetic Bragg peaks, implying the absence of 
ordinary magnetic long-range order. As for the static 
spin structure factor 5(k), the three-sublattice nematic 
order shows a sharp peak at the corners of the re- 
duced Brillouin zone that are inside the original Brillouin 
zone of the triangular Bravais lattice, and the spin-spin 
correlation function shows a power-law decay in space 
(S(r)S(O)) — l/|r| 3 . We should emphasize that the peak 
in S(k) does not diverge but has only a kink singularity, 



<S(k) ~ const — A|k — ko|. In neutron experiments, the 
spin correlation length was determined as the inverse of 
the peak width, and therefore that is consistent with our 
result. There remain two points that should be further 
investigated in future studies. The first point is a detailed 
structure of S(k). The peak position is different between 
the neutron data and the theoretical calculation shown 
before. We believe that it is possible to reproduce similar 
S(k) by tailoring the starting Hamiltonian by including 
longer-range interactions, like in the case of incommen- 
surate helical spin order. We should emphasize that the 
basic features presented above will not change aside from 
detailed structure in S(k). The second point is the energy 
dependence. We also need more quantitative analysis on 
the physical properties at finite temperatures and also 
under finite magnetic fields, but we believe that the sce- 
nario presented in this paper will help the understanding 
of the intriguing properties of NiGa2S4. 
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